An ab-initio many-body method for electronic structure calculations of solids. 

I. Description of the method 
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We propose a new, alternative method for ab-initio calculations of the electronic structure of solids, 
which has been specifically adapted to treat many-body effects in a more rigorous way than many 
existing ab-initio methods. We start from a standard band-structure calculation for an effective 
one-particle Hamiltonian approximately describing the material of interest. This yields a suitable 
set of one-particle basis functions, from which well localized Wannier functions can be constructed 
using a method proposed by Marzari and Vanderbilt. Within this (minimal) basis of localized 
Wannier functions the matrix elements of the non-interacting (one-particle) Hamiltonian as well 
as the Coulomb matrix elements can be calculated. The result is a many-body Hamiltonian in 
second quantization with parameters determined from first principles calculations for the material 
of interest. The Hamiltonian is in the form of a multi-band Hamiltonian in second quantization 
(a kind of extended, multi-band Hubbard model) such that all the standard many-body methods 
can be applied. We explicitly show how this approach can be solved in the simplest many-body 
approximation, the mean-field Hartree-Fock approximation (HFA), which takes into account exact 
exchange and corrects for self-interaction effects. 

PACS numbers: 71.10.Fd, 71.15.AP, 71.15.Mb, 71.20.Be, 71.45.Gm, 75.10.Lp 



I. INTRODUCTION 

Most existing ab-initio (first-principles) methods for 
the numerical calculation of the electronic properties of 
solids are based on density functional theory (DFT) 1 , 
which in principle is exact and properly takes into ac- 
count many-body effects involving the Coulomb interac- 
tion between the electrons; for an overview on the present 
status of DFT we refer to the booka 2 ^. But, in gen- 
eral, the functional dependence of the kinetic energy and 
the exchange and correlation part of the Coulomb (in- 
teraction) energy on the electron density are not known 
explicitly, and hence additional approximations and as- 
sumptions are necessary. A well established additional 
approximation is the local density approximation (LDA) 4 
(or local spin-density approximation, LSDA, for magnetic 
systems) , which assumes that the electron density of the 
interacting system is the same as that of an effective 
non-interacting system and that the exchange-correlation 
potential depends only on the electronic density locally. 
Even then, the functional dependence of the exchange- 
correlation energy on the density is not known in general, 
and it is usually necessary to make an ansatz for the 
exchange-correlation functional, which is based on the 
homogeneous electron gas. The LDA goes beyond the 
simplest electron-gas approximation, the Hartree-Fock 
approximation (HFA), in that correlation energy (i.e. the 
part of the interaction energy beyond HFA) is explic- 
itly taken into account. More recent generalized gradi- 
ent approximation (GGA) calculations also make some 
improvements in the correlation energy that slightly im- 
prove on the homogeneous electron gas approximation. 
On the other hand, the exact HFA exchange potential 
is non-local, i.e., it depends on the electron wave func- 



tions and density at all other positions, an effect which 
the local LDA exchange potential misses. Also, LDA 
misses some cancellation between the exchange and cor- 
relation parts of the total energy (self-interaction correc- 
tions). However, in practice LDA-treatments are sim- 
pler than HFA-calculations (especially for metals), be- 
cause local exchange is easier to treat than non-local 
exchange, and are usually in better agreement with ex- 
periment. Therefore, DFT-(LDA-like) treatments have 
been far more common than HFA during the past few 
decades, even in quantum chemistry (with a long tradi- 
tion of methods based on HFA). 

Ab-initio DFT calculations based on the L(S)DA have 
been very successful for many materials and ground-state 
properties such as crystal structure, ground state and 
ionization energy, lattice constant, bulk modulus, crys- 
tal anharmonicity^, magnetic moments, and some photo 
emission spectra. However, there are also important lim- 
itations. For example, LDA predicts a band gap for semi- 
conductors that is almost a factor of two too small. In 
addition, for many strongly correlated (narrow energy 
band) systems such as high-temperature superconduc- 
tors, heavy fermion materials, transition-metal oxides, 
and 3d itinerant magnets, the LDA is usually not suf- 
ficient for an accurate description (predicting metallic 
rather than semiconducting behavior, failing to predict 
quasi-atomic- like sattelites, etc). 

Therefore, it is justified and important to look for 
new, better ab-initio methods and improvements that go 
beyond L(S)DA. In fact, there have already been sev- 
eral attempts to improve LDA, such as gradient correc- 
tions, non-local density schemes, self-interaction correc- 
tions, and the GW approximation. Gradient corrections 6 
approximately account for the fact that the electron den- 
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sity is not constant but r-dependent in an inhomoge- 
neous electron gas; estimates on the influence of pertur- 
bations of the homogeneity on the ground state energy 
of the homogeneous electron gas are then used to con- 
struct a new functional for the exchange-correlation po- 
tential containing Vn(r) terms. The non-local density 
schemes go beyond LDA by considering that the exact 
exchange-correlation potential V xc (r) cannot depend only 
on the density n(r) at the same position r but should 
depend also on the electron density at all other posi- 
tions n(r'). Usually the new ansatz for the functional 
of the exchange-correlation energy contains the pair cor- 
relation function or the interaction of the electrons with 
the exchange-correlation hole Ai. The recently developed 
exact exchange (EXX) formalism- cancels the spurious 
(unphysical) electronic self-interaction present in LDA 
and gradient corrected exchange functionals. A stan- 
dard method for ab-initio calculations of excited states 
is the GW-approximation (GWA) 9,10 . Denoting the one- 
particle Green function by G and the screened interac- 
tion by W, the GWA is an approximation for the elec- 
tronic self energy £ « GW , which is correct in linear or- 
der in W and can diagramatically be represented by the 
lowest-order exchange (Fock) diagram. The one-particle 
Green function G is usually obtained for the effective 
one-particle LDA Hamiltonian. 

Recently there have been several attempts to 
combine ab-initio LDA calculations with many-body 
approximations, 11 ! 12 ! 13 ! 14 ! 15 ! 16 ! 17 ! 18 ! 19 ! 20 All of these re- 
cent developments add local, screened Coulomb (Hub- 
bard) energies U between localized orbitals to the one- 
particle part of the Hamiltonian obtained from an ab- 
initio LDA band-structure calculation, but differ in how 
they handle the correlation part. In the earliest attempts, 
the LDA+U methodii used essentially a static mean- 
field-like (or Hubbard-I-like) approximation for the cor- 
relation. The simplest approximation beyond Hartree- 
Fock, second-order perturbation theory (SOPT) in U, 
was usediSii 3 ^ 6 ^ 9 " to study the electronic properties of 
3d-systems (like Fe and Ni) and heavy fermion sys- 
tems (like UPt 3 ). The LDA++ approach 15 ! 17 ; 18 has a 
similar strategy, but uses other many-body approxima- 
tions to treat the correlation problem, namely, either the 
fluctuation-exchange approximation (FLEX) or the dy- 
namical mean-field theorySi (DMFT); for a review on 
this recently very successful LDA+DMFT 20 approach see 
Ref. [23. These approaches, including the LDA+U, have 
in common that they have to introduce a Hubbard U as 
an additional parameter and hence are not really first- 
principles (ab-initio) treatments. Although they use an 
LDA ab-initio method to obtain a realistic band struc- 
ture, i.e., single-particle properties, Coulomb matrix el- 
ements for any particular material are not known, and 
the Hubbard U remains an adjustable parameter. In ad- 
dition, since some correlations are included in LDA as 
well as by the Hubbard U, it is unclear how to separate 
the two effects and double-counting of correlation may 
be included in these approximations. 



In this paper we suggest a different, alternative ap- 
proach and propose a new type of ab-initio treatment. 
Because we want to avoid DFT-LDA due to its problems 
with double counting and self-interactions, we instead 
suggest starting from a realistic Hamiltonian in second 
quantization and directly applying many-body methods 
to this Hamiltonian. The underlying many-body the- 
ory in DFT-LDA involves the homogeneous electron gas, 
from which a functional for the exchange-correlation po- 
tential is derived. But we want to apply the many-body 
theory to a Hamiltonian tailored to the material of in- 
terest, i.e., study a problem of interacting electrons on 
a lattice and therefore an inhomogeneous electron gas. 
To get the realistic, microscopic Hamiltonian in second 
quantization the one-particle and Coulomb matrix ele- 
ments have to be calculated. For that purpose one needs 
a suitable one-particle basis. For a solid it is reasonable 
that either Bloch or Wannier functions should form a 
suitable one-particle basis set. We suggest starting from 
a set of maximally localized first-principles Wannier func- 
tions, because the good localization of these wave func- 
tion should make the one-particle (tight-binding) and 
two-particle (Coulomb) matrix elements important only 
for a finite number of near-neighbor shells (on-site, near- 
est, next nearest neighbors, etc.) such that only a few 
matrix elements have to be actually calculated. It is 
possible to determine these Wannier states for an ef- 
fective (auxiliary) one-particle Hamiltonian. For exam- 
ple, one could use the true non-interacting Hamiltonian, 
the Hartree Hamiltonian, or even the LDA Hamiltonian. 
In principle, the final results should not depend on the 
choice of this auxiliary one-particle Hamiltonian, because 
any set of Wannier functions should span the whole one- 
particle Hilbcrt space. But, in practice, one typically 
works only within a restricted (finite dimensional) sub- 
space. Therefore, the auxiliary one-particle Hamiltonian 
and the resulting one-particle basis need to be chosen 
in such a way that the physically most important sub- 
space is spanned. In practice, we have found the Hartree 
Hamiltonian to be a suitable choice for the auxiliary, ef- 
fective one-particle Hamiltonian. 

Our suggested procedure is schematically outlined as 
follows: First, a traditional band-structure calculation is 
used for the auxiliary one-particle Hamiltonian to calcu- 
late the eigenstates in the form of Bloch functions. Wan- 
nier functions are closely related to these Bloch functions 
via a unitary transformation, and thus span the same 
one-particle space as the Bloch functions. However, since 
the phases of the Bloch functions are arbitrary, Wannier 
functions are not unique. Their non-uniqueness (gauge 
freedom) is used to construct "maximally localized Wan- 
nier functions" using a method proposed by Marzari and 
Vanderbilt^ 3 .. A proper localization of the Wannier func- 
tions is important, because only then do the standard 
assumptions that are frequently used in model treat- 
ments hold, e.g., that both one-particle (tight-binding) 
and two-particle (Coulomb) matrix elements are impor- 
tant only on-site and for a few neighbor shells. Next, 
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in order to describe the physical Hamiltonian in second 
quantization, one has to evaluate a restricted number of 
hopping and Coulomb matrix elements within the ba- 
sis of the well localized Wannier functions. Making use 
of the special representation of the wave functions in- 
herent to the specific band-structure method used (for 
instance in terms of spherical harmonics, plane waves or 
Gaussian orbitals) can simplify the actual computation 
of these matrix elements; for the Coulomb matrix ele- 
ments we propose, in particular, an algorithm based on 
fast Fourier transformation, which works independent of 
the special band-structure method. The result of this 
procedure is a multi-band, second-quantized Hamilto- 
nian in a Wannier representation (a kind of extended, 
multi-band Hubbard model) describing interacting elec- 
trons on a lattice. Unlike model studies of correlated 
lattice electrons, the matrix elements of this Hamilto- 
nian are not free adjustable parameters, but are obtained 
from first principles for the given material. This Hamil- 
tonian should be evaluated by an appropriate many- 
body methods. The simplest approximation is the mean- 
field Hartree-Fock approximation (HFA), which exactly 
treats exchange and excludes self-interaction effects. Un- 
like first-quantization treatments of HFA, which require 
the solution of a Schrodinger equation with a non-local 
potential, in second-quantization only the selfconsistent 
evaluation of expectation values is required (once an ap- 
propriate one-particle basis has been chosen). But, of 
course, HFA is certainly not sufficient for realistic materi- 
als properties; higher order correlations have to be taken 
into account, for example, perturbationally by resumma- 
tions of Feynman diagrams (RPA-bubble diagrams lead- 
ing to screening, ladder diagrams, etc.) or employing 
recently developed non-perturbational methods like dy- 
namical mean-field theory (DMFT) (if the basic assump- 
tions, for instance concerning the locality of the most 
important interaction terms, hold for the specific many- 
body Hamiltonian). 

This paper is organized as follows: In Section ^ we 
shortly repeat some basic notations, give the Hamilto- 
nian in first and second quantization, describe DFT and 
LDA, Wannier- and Bloch-states as possible basis states 
spanning the one-particle Hilbert space and the necessary 
truncation. Section inTl describes how a suitable restricted 
basis set is obtained by solving the Schrodinger equation 
for an auxiliary one-particle Hamiltonian, which should 
be chosen appropriate for the material of interest and in- 
corporating information about the lattice structure, lat- 
tice constant, total electron number, etc. In Section llVl 
we discuss how another equivalent basis set of maximally 
localized Wannier functions can be obtained from this 
one-particle basis. Within this Wannier basis ab-initio 
one-particle (tight-binding) and two-particle (Coulomb) 
matrix elements can be calculated, as outlined in sec- 
tion The resulting many-body Hamiltonian in second 
quantization with ab-initio parameters is then ready to 
be studied by suitable many-body approximations, which 
is the subject of Section IVll We then discuss and com- 



pare the advantages and disadvantages of our proposed 
ab-initio method relative to other ab-initio methods in 
Section lYlII The following paper (II) contains an initial 
simple application (an unscreened HFA) of this approach 
to the 3d transition metals Fe, Co, Ni and Cu. 



II. HAMILTON OPERATOR IN FIRST AND 
SECOND QUANTIZATION 

A system of N e interacting (non-relativistic) electrons 
can be described by the Hamiltonian 



N„ 



N s 2 

i—l i—1 i>j J 



(i) 

The first part T is the kinetic energy of the electrons. 
The V(r) describes the external one-particle potential, 
which for molecules and solids is solely caused by the 
positive nuclei or ions and-within the adiabatic or Born- 
Oppenheimer approximation-can be regarded as a static 
potential. For solids, i.e., for electrons in a crystal, V(r) 
is a periodic potential. The third part W is the Coulomb 
interaction between the electrons. 

For the interacting system it is impossible to exactly 
solve the full Schrodinger equation and to calculate the 
many-body eigenenergies (or even the ground state en- 
ergy) and the totally antisymmetric many-particle wave 
functions, which can be represented as a linear combina- 
tion of Slater determinants (a configuration-interaction 
or CI approach to correlation). This is practical only 
for a very small number of electrons. The basic idea of 
DFT is to avoid determining the antisymmetric (ground 
state) wave function and instead to focus on the (ground 
state) density n(r). This requires knowing the functional 
dependence of the kinetic energy T, the one-particle po- 
tential V, and the interaction energy W on the den- 
sity n. From this the ground state density can be de- 
termined from a variational principle. But in general 
this DFT-concept can only be approximately applied, 
since the functional dependence of the kinetic energy 
and the exchange-correlation energy on the density is un- 
known. Therefore additional assumptions and approxi- 
mations are necessary, such as the LDA where the den- 
sity n(r) is the same as a non-interacting electron system 
in an effective one-particle potential and the exchange- 
correlation energy depends only locally on the electron 
density. From a fundamental point of view these addi- 
tional assumptions are hard to justify and uncontrolled, 
and an estimate of the error or systematic corrections are 
very difficult to obtain. 

In this paper we choose yet another approach that em- 
ploys the (very elegant) formalism of "second quantiza- 
tion" , which automatically accounts for the antisymme- 
try through the fermion anticommutation relations. In 
second quantization the full many-body Hamiltonian can 
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be written as: 
H ' = ^2 t ij c h c J° + \ W ljM c\ a c] al c k<jl ci (J (2) 

i,j,<T i,j,k,l,a,tr' 

Here i,j,k,l denote a complete set of one-particle or- 
bital quantum numbers and c, a 1 are the spin quantum 
numbers. The states \i) and the corresponding wave func- 
tions (fii(r) = (r\i) form a basis of the one-particle Hilbert 
space. The matrix elements in Eq. |21are defined by: 

= J d 3 r^(r)(-^V 2 + V(r)) Vi (r) 

w ijM = m\T^\m) (3) 

= J d\ J d 3 rV:W^(r')^^(r')w(r) 

Of course these matrix elements depend on the one- 
particle basis {\i}} that is chosen. But because of the 
completeness relation the physical results should, in prin- 
ciple, not depend on the choice of the one-particle basis. 
Because of the lattice periodicity an obvious choice for a 
one-particle basis is a Bloch basis {|nk)}; then the orbital 
one-particle quantum numbers n, k are the band index n 
and the wave number k (within the first Brillouin zone). 
The Bloch wave functions can be written as a product of 
a plane wave with a lattice periodic Bloch factor 

V>nk(r) = (r|nk) = e ikr w„ k (r) (4) 
Unk(r) = M„k(r + R) 

where R denotes a lattice vector. In practice one can 
work only on a truncated, finite-diemnsional one-particle 
Hilbert space. Here the truncation consists of including 
only a finite number of bands and a set of k-values from 
a discrete mesh in k-space. A Bloch basis can be ob- 
tained by applying a traditional band-structure method 
to solve the Schrodinger equation for an effective one- 
particle Hamiltonian with a periodic potential. But, 
because the Bloch states are delocalized, a very large 
number of Coulomb matrix elements (depending on four 
quantum numbers) between all possible k-states would 
have to be calculated. Therefore, it seems that a more 
appropriate basis would be to use well localized wave 
functions. Although the one-particle part is not diago- 
nal within such a localized basis, it is expected that a 
short-range tight-binding assumption will hold, i.e., that 
the on-site and the intersite matrix elements for only a 
few neighbor shells are sufficient. It is well known from 
elementary solid state theory that Wannier states pro- 
vide an alternative orthonormal basis set, which spans 
the same space as the Bloch states. The Wannier states 
are related to the Bloch states by the unitary transfor- 
mations: 

w Rn (r) = (r|Rn) = l^e- lkR V„ k (r) (5) 

k 



l^nk) = ]Te* kR |Rn) 

R 

Now our strategy is the following: 

• Perform a traditional band-structure calculation 
for an effective one-particle Hamiltonian H c ff with 
lattice periodicity to obtain a Bloch basis of the 
Hilbert space. Only a finite number of band indices 
will be considered and the calculations will be done 
for a discretized, finite mesh in k-space, i.e., we will 
work only on a reduced, truncated Hilbert space. 

• Determine well-localized Wannier functions span- 
ning the same (truncated) Hilbert space as the 
Bloch basis from the canonical transformation (J5J 
described above. 

Of course, the important energy bands (and correspond- 
ing band indices) are those that determine the electronic 
properties of the system, i.e., the bands near to the Fermi 
level. Because the Hilbert space is truncated, we do no 
longer work with a complete basis set. Hence, it is im- 
portant to start from Bloch functions obtained from a 
band-structure calculation for a well chosen effective one- 
particle Hamiltonian. Possible and suitable choices for 
the effective (auxiliary) one-particle Hamiltonian will be 
discussed in the next Section lTTTl Furthermore, the Wan- 
nier functions obtained according to JSJ are not unique, 
because the phases of the Bloch functions can be cho- 
sen arbitrarily. In fact, for a given set of randomly cho- 
sen phase factors the Wannier functions may not be lo- 
calized at all. On the other hand, one may use this 
gauge freedom to construct optimally localized Wannier 
functions 23 , as will be described in Section llVl 



III. BAND STRUCTURE CALCULATION FOR 
THE AUXILIARY ONE-PARTICLE PROBLEM 

Our Bloch state basis should be suitably chosen for 
the material of interest, and should be obtained from a 
standard band-structure calculation of an effective one- 
particle Hamiltonian with a lattice periodic potential: 

V eG (r) = V cS {r + R) (6) 

for any lattice vector R. The corresponding one-particle 
Schrodinger equation is 

2 

(|-+Kff(r))Vnk(r) = £„(k)V„k(r) . (7) 

Different choices for the effective periodic potential are 
possible. The simplest choice would be the bare one- 
particle potential V(r). However, since the Coulomb 
interaction is strong, in general, the resulting eigenen- 
ergies of this non-interacting Hamiltonian will be much 
lower than the relevant energies of the interacting sys- 
tem. For example, in the 3d-systems like Ni or Cu the 
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3d- and 4s-states form bands close to the Fermi energy. 
However, without any Coulomb repulsion the 3d-states 
become very strongly bound atomic-like (core) states, 
which would be pushed well below the Fermi energy, and 
therefore the corresponding Bloch cigcnfunctions are not 
a good starting point to describe the electronic bands 
close to the Fermi level. Because the Hilbert space is 
truncated, it is extremely important to start from a band 
Hamiltonian T + V e & that gives eigenfunctions as close as 
possible to those which are expected to form the relevant 
many-body states of the interacting system. The bare 
one-particle potential is consequently a bad choice. 

Another possible choice for the effective periodic one- 
particle Hamiltonian is the Hartree Hamiltonian. Then 
the Bloch basis is obtained by solving the one-particle 
Schrodinger equation 



't- + V ( r ) + Vk(r)W(r) = e„(k)^„ k (r) (8) 
vim / 

where the Hartree potential is given by 



(9) 



and the particle density 

p(r)=£/(ei(k))llMr)| 



(10) 



Another possible choice for the effective (auxiliary) one- 
particle potential would be the LDA-potential. Because, 
in this context, LDA would only be used to construct 
basis functions, there would be no problem with double 
counting interaction terms. 

After choosing the appropriate effective one-particle 
potential one has to employ a numerical method to solve 
the Schrodinger equation J7J for the periodic potential. 
For this established and well known band-structure meth- 
ods can be used, such as the "linearized muffin-tin or- 
bital" (LMTO) method 2 ^, the (full potential) "lin- 
earized augmented plane wave" (LAPW) method (us- 
ing the WIEN2k computer code^I) , or the "full-potential 
local-orbital" (FPLO) 2 ^ minimum-basis method. 



IV. DETERMINATION OF MAXIMALLY 
LOCALIZED WANNIER FUNCTIONS 

The band-structure calculation diagonalizes the aux- 
iliary one-particle Hamiltonian to determine the Bloch 
states V'nk(r) = e lkr u n k(r) labeled by band indices {n}. 
After truncating this one-particle Hilbert space by con- 
sidering only a finite number J of band indices, unitary 
transformations are possible for each fixed k that lead to 
a new basis set 



Ika 



(13) 



where f(E) 



1 



3 0(E-fj.) _|_ i 



(11) 



is the Fermi function (step function for the ground state 
T = or (3 = oo). Since the occupied states \ipik) 
determine the particle density and Hartree potential, 
Eq. [S]has to be solved selfconsistently, which is usually 
done by iteration procedure. The advantage of includ- 
ing the Hartree potential into the effective potential, i.e. 
T4fr(r) = V(r) + Vk(r), is that the approximate effects 
of the Coulomb interaction are included (in the mean- 
field approximation). Therefore, the eigenenergies (en- 
ergy bands) will be about the right magnitude and the 
resulting basis functions can be expected to be more suit- 
able in the energy regime around the Fermi level. 

Since the only purpose in solving the effective one- 
particle Schrodinger equation J2J is the construction of 
a suitable basis set of Bloch functions, we will not make 
use of the eigenenergies e„(k) obtained in J7J) or give 
these solutions any physical interpretation. For this rea- 
son one could also use a variety of different artificial effec- 
tive potentials. For example, we have studied a "weighted 
Hartree potential" 



V bB (t) = V(t)+xV r {t) 



(12) 



with an additional parameter x, which can be varied 
so that the choice of the basis functions minimizes the 
resulting Hartree- or Hartree-Fock ground state energy. 



which are still Bloch functions but with other band in- 
dices (within which the auxiliary one-particle Hamilto- 
nian is, in general, no longer diagonal). One can use this 
gauge freedom to construct a Bloch basis (still of the 
same restricted one-particle Hilbert space) for which the 
corresponding Wannier functions determined according 
to JSJ are maximally localized. 

A suitable measure for the localization of the Wannier 
function is the spread functional 



n = 



[(r 2 )n (r)l] 



(14) 



where the notation (A) n = (On\A\On) for any operator 
A has been used. Among all the possible gauges, we will 
use that gauge which minimizes this spread functional Q. 

A method for minimizing Eq. (|14H has been developed 
by Marzari and Vanderbilt 23 -. This method has already 
found widespread applications recently 29,30 . It starts 
from a decomposition into invariant, diagonal, and off- 
diagonal contributions: 



E 



2 }„-^l(Rm|r|0n) 



Rn 



(15) 



£ \(*-rn\r\0n) 
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The first term is "gauge- invariant" , i.e., independent of 
the choice of unitary transformations among the bands. 
The second term can be decomposed into band-off- 
diagonal and band-diagonal components: 



n = 



^2^2\{Rm\r\On)\ 

m^n R, 



■EE 



(Rrt|r|On) 



OD 



Or 



(16) 

In practice one needs the unitary matrices U^ m on a dis- 
cretizcd k-mesh, which can be chosen to be simple cubic. 
If b denotes the vectors connecting each k-point to its 
nearest neighbors, one can define 



M mn = (WmkK,k+b) 



(17) 



where u„k are the Bloch factors appearing in ijlfl. Using 
(|17|l one can express the expectation values of r and r 2 
asSi: 



^E^ b 



Im In M. 



k.b 



W b 



k.b 



w (18) 



( r2 )™ = ^E^{[l-|^ b r] + PmlnM^ b ] 2 } 

k.b 

These expressions are not unique. This non-uniqueness 
arises from the fact that it is only required that these ex- 
pectation values become exact in the limit of a dense 
mesh. The present choice guarantees that under the 
gauge transformation u„k) — > e~ lkR |it n k) (correspond- 
ing to a translation by a lattice vector) the following 
physical properties are fulfilled: 



(r)r 



<r)„ + R 



(19) 



2(r)„R + i? 2 



Then one finds for O = Ot + O d + Or 



OD 



Sir 



k.b \ mn / 

^E^E «i 2 

k,b m^n 

4 E w " E (- Im ln M - - b < r >«) : 



(20) 



k.b 



The Marzari-Vanderbilt algorithm 2 ^ now consists in an 
iterated application of a (in principle infinitesimal, in 
practice numerically discrete) canonical transformation 
chosen so that the gradient of O is negative, i.e., the 
spread functional and thus the derealization of the Wan- 
nier function decreases. In each iteration cycle Oj should 
not change whereas OD + Od should decrease. 

In practice the iteration cycle can schematically be 
written down aa2&: 

C/W k = [/( W -D k exp [AW k (M( JV - 1 ) k ' b )l (21) 



M (AT)k,b = ^AQk f M (JV-l)k,b jj(N)k+b ^2) 

with the initialization: 

(23) 



u mn ~ u mn 



M£l k ' h = (u mk \u nM+b ) 



(24) 



Here AIU k (Af k ' b ) is an anti-hermitian matrix (so that 
exp [AIU k ] is unitary) which according to Ref. |23 is 
explicitly given by the following equations: 



AIU k = - w b (A(R* h ) - S(T kb )) 



(25) 



where A{B) = (B - B^)/2,S(B) = (B + B^)/2i are 
the symmetrizing and antisymmetrizing operations for 
the operators or matrices B, w — ^^Wb, ol 6 [0,1] is a 
numerical parameter determining the discrete, finite step 
size, and the J x J-matrices i? kb and T kb are explicitly 
given by: 



R 



kb 

mn 



M kb M kb * 

mn nn 



(26) 



-kb 



M 



kb 



M kb 

nn 



(lmlnM kb + b(r)„) 



The whole algorithm (|21J) can be considered as an itera- 
tion scheme to construct the unitary transformation (7 k . 
Once this iteration has converged one can perform the 
transformation to the new Bloch functions according to 
(|13fl and then calculate the optimally localized Wannicr 
functions according to the definition JSJ. 

In practice it is useful to prepare the Bloch orbitals to 
make the starting Wannier functions somewhat localized. 
This has two advantages: (i) the minimization proce- 
dure converges faster, and (ii) this helps to avoid getting 
trapped in local minima. For that purpose we make a 
rather trivial "gauge transformation" to each Bloch func- 
tion obtained by the band structure calculation, namely 
a multiplication with a phase factor according to 



^nk(r) 



exp [ - i Im ln V>nk( r n)] ipnk 



(r) ■ (27) 



This gauge transformation has the property that 
Im ln^ , „k(r n ) transforms to zero. So at the point r n , all 
the Bloch functions will have the same phase and (r„|0n) 
will take a large value. To make the method work well, 
one should choose r„ where the Wannier functions are 
expected to be reasonably large and one can choose r„ 
individually for each band n. 



V. ONE-PARTICLE AND COULOMB MATRIX 
ELEMENTS 

After the maximally localized Wannier functions 
%„(r) have been determined, the next task is to cal- 
culate the one-particle and Coulomb matrix elements of 
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the Hamiltonian. The one-particle matrix elements in 
Wannier representation are explicitly given by: 



d 3 rw$(r) 



2m 



V 2 + V(r) w 2 (r) (28) 



where V(r) is the lattice periodic one-particle potential 
for the non-interacting electrons. Here, and in the fol- 
lowing, we use the abbreviated notation 1 to mean Rini 
and 2 to mean for R2^2, etc. Note that the one-particle 
part of the Hamiltonian cannot be expected to be band 
diagonal because of the unitary transformation to an- 
other set of the band indices within which the Wannier 
functions are maximally localized. But by construction 
and definition the Wannier functions obey the following 
translational invariance property: 



ifRn(r) = w „(r - R) 



(29) 



where is the arbitrarily chosen zero lattice vector. 
Thus, the one-particle matrix elements depend only on 
the relative distance R = Ri — R2 : 



(Rnl 



2m 



V 2 + V(r) \0m) (30) 



Because of the strong localization of the Wannier func- 
tions one can safely expect that these one-particle (hop- 
ping) matrix elements decrease with increasing distance 
|R| and have to be evaluated explicitly only for a limited 
number of R, in particular on-site, i.e. for 0, and for 
a few neighbor shells R. In practice it turned out that 
for properly localized Wannier functions (optimized ac- 
cording to the Marzari-Vanderbilt algorithm) the explicit 
evaluation of one-particle matrix elements up to the 5th 
neighbor shell is sufficient, whereas for non-optimally- 
localized Wannier functions (made only "somewhat lo- 
calized" according to the prescription l|27|) , for instance) 
up to 30 neighbors have to be taken into account i^i De- 
pending on the specific band-structure method used for 
the effective (auxiliary) one-particle problem, the eigen- 
functions (Bloch functions) are usually represented as lin- 
ear combinations of certain elementary functions, for in- 
stance plane waves, Gaussians, or spherical harmonics. 
Thus, also the evaluation of the three-dimensional inte- 
gral in (|28f) can usually be traced back to known integrals 
over these elementary functions, and usually at most a 
one-dimensional integral remains to be calculated explic- 
itly numerically. 

The Coulomb matrix elements are given by: 

^12,34= [ d 3 vd 3 v' <(r) W *(r') — 6 — w 3 {v') w A (r) 
J |r-r'| 

(31) 

Let us first take a brief look at general properties of the 
matrix elements which are useful for minimizing com- 
puting time and memory storage. From i|31|) and l|29() . it 
follows that 



W 



12,34 



W, 



(R 1 -R 4 ?i 1 )(R 2 -R 4 ri2),(R.3-R-4™3)(On4) 



(32) 



That is, we may always translate the lattice site indices 
in a way that R 4 — > 0. Moreover, since r and r' in 
(ffi"Jl can be interchanged, we have ^12,34 = W2i,43 and 
since W is Hermitian we have Wi2,34 = WJ 3 21 . For the 



practical evaluation of these 6-fold integrals a fast Fourier 
transformation (FFT) algorithm turned out to be very 
efficient and independent of the specific representation of 
the Wannier functions. Using 



d 3 q 



e lqr 2?r 2 



one obtains: 

^12,34 

where 

/i(q) = 



7^2 /^ 3 q \ A(q) M-q) 

Itt z J q z 



(33) 



(34) 



! re iq X ini (r) %M (r) (35) 



i2(q) = |rf ; 



r e 



y R 2 "2 



(r) % 3 n 3 ( r ) ( 36 ) 



These functions are just the Fourier transforms of a prod- 
uct of Wannier functions. They can be calculated very 
efficiently by evaluating the Wannier functions on a cu- 
bic mesh in r-space with some finite spacing Ax and then 
applying a standard (numerical) FFT-algorithm. The re- 
sult of this Fourier transformation is /i.2(q) on a cubic 
mesh in q-space with some Aq. As F(q) = /i(q) fi{— q) 
is a smooth function at q = 0, a divergence arises in the 
integrand from q~ 2 . In order to treat this divergence, we 
split the integral by subtracting and adding F(0): 

q J q J q 

(37) 

The first non-vanishing term of a polynomial expansion 
of the numerator F(q) — F(0) is of order q 2 . Hence, 
the polynomial expansion of the integrand starts with a 
constant and the divergence is avoided. 

All integrals are over a cube with length 2p — NAq. 
The first integral in (|37|l is evaluated by transforming the 
integral into a sum over little cubes with volume (Aq) 3 . 
At q = 0, the value of the integrand is calculated via the 
second derivative of F(q) at q — numerically 

The remaining integral in l|37|) is simply a constant 
given by: 



dq x 



-p 



dq % 



P 



• I dq * 7a 

dq x / dq y / dq z — 

1 J-i J -1 



V ■ C (38) 



with C = 15.34825, which we have evaluated numerically. 

As argued already for the one-particle matrix elements, 
the Coulomb matrix elements need explicitly be evalu- 
ated only on-site, i.e. for Ri = R2 = R3 = 0, and for 
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at most only a few neighbor shells, probably only nearest 
neighbors, because of the good localization of the Wan- 
nier functions used. Here again the application of the 
Marzari-Vanderbilt algorithm turns out to be important 
in order to reduce the number of explicitly necessary eval- 
uations of Coulomb matrix elements, which are (compu- 
tational) time consuming. 

VI. MANY-BODY TREATMENT OF SECOND 
QUANTIZED HAMILTONIAN WITH AB-INITIO 
PARAMETERS 

After the evaluation of the one-particle and two- 
particle matrix elements we have the Hamiltonian in sec- 
ond quantization, i.e., in the form (J2J in a Wannier repre- 
sentation with all parameters determined from first prin- 
ciples. The only approximations made so far are the re- 
striction to a finite number J of band indices for a specific 
basis set and the approximations inherent to the band- 
structure calculation method used (for example, muffin- 
tin potentials, linearization approximations, discretiza- 
tion in k- and real r-space, etc.). The Hamiltonian is in 
the form of a multi-( J-)band extended Hubbard model. 
A one-band description is usually not sufficient even for 
the simplest solids. Furthermore, the hopping matrix el- 
ements are in general not restricted to nearest neighbors 
(but taken into account up to the 5th nearest neighbors) 
and direct and exchange Coulomb matrix elements on- 
site, inter-band and intra-band and, if necessary, also 
inter-site (but again only up to a few neighbor shells) 
are taken into account. 

Now, in principle, one can apply any many-body 
method that has been developed for interacting electron 
systems on a lattice. Many of the standard methods rely 
on systematic perturbation theory with respect to the 
Coulomb interaction, and the terms can be represented 
by Feynman diagrams^ 2 .. The lowest-order selfconsistent 
approximation within this diagrammatic approach is the 
Hartree-Fock approximation (HFA), which can also be 
derived from a variational principle. Within the HFA 
the selfenergy is explicitly given by: 

vHF v^Hart , v^Fock /QO\ 

^12,(7 — ^12 + ^12,0- l dy J 

= Yl 1^13,42 - <WW31,42] (4 ct ,C 4ct ') 
34ct' 

Here the expectation values (or the density matrix) 

A\ 2 = (c\ a c 2a ) (40) 

have to be determined selfconsistently for the Hartree- 
Fock Hamiltonian: 

tfHF = J2 (* 12 + <L<** ( 41 ) 

12a 

These expectation values are easily calculated by going 
back to the Bloch states within which the effective one- 
particle Hamiltonian (|41|l is diagonal with eigenenergies 



E° n (k) . In this Bloch representation the expectation val- 
ues for zero temperature (i.e., within the ground state) 
are given by 

(cj^cvkv) = S nn ,6 kk ,S aa ,6{E F - K(k)) . (42) 
One finds 

A i2 = Jf E e - lk(Rl ~ R2) £ U* ni U** n2 9(E F ££(k)) 

k m 

(43) 

where (similar as in the U^ n denote the unitary 

transformation between the basis with band indices lead- 
ing to maximally localized Wannier functions and the 
basis within which the effective one-particle Hamiltonian 
H41JI is diagonal. For the total energy in the Hartree-Fock 
approximation one finds: 

Euf = = £ (*u + \^ + ^Sg?) ^12 (44) 

This follows formally from the original many-body 
Hamiltonian J2J using the decoupling 

(45) 

which is valid according to Wick's theorem, if the ex- 
pectation values are calculated with respect to an effec- 
tive one-particle Hamiltonian, which the Hartree-Fock- 
Hamiltonian 1)41(1 is. 

Of course the HFA is only the simplest and basic many- 
body approximation which can be applied (mean-field). 
It is correct only up to linear order in the Coulomb in- 
teraction. The main advantage of the diagrammatic ap- 
proach to the second-quantized many-body Hamiltonian 
P|) is that one can formulate systematic improvements 
and corrections to the basic approximations. One possi- 
bility to go beyond HFA is to take into account all skele- 
ton diagrams up to second order in the Coulomb inter- 
actions. Probably better and more reliable many-body 
approximations are obtained by resummations of an in- 
finite series of diagrams of a certain class. For instance 
the resummation of all bubble diagrams corresponds es- 
sentially to the random phase approximation (RPA) and 
means that in the Fock diagram the bare interaction has 
to be replaced by an effective screened interaction. 

Of course the application of many-body theory is not 
restricted to the standard perturbational methods in 
terms of Feynman diagrams. During the last decade 
much progress in many-body theory has been made due 
to the development and successful application of non- 
perturbative methods, such as the DMFT—. In its orig- 
inal version the DMFT requires a local (k-independcnt) 
selfenergy, which is usually justified for 3-dimcnsional 
systems. But even if this assumption should not be jus- 
tified, the application of non-local (cluster) extensions of 
DMFT is possible. In any case, for realistic materials a 
multi-band system has to be studied; therefore a mapping 
on a multi-level single-impurity problem and new types 



9 



of DMFT-selfconsistency cycles are probably required. 
Another possible non-perturbative many-body approach 
to our Hamiltonian J2J) with ab-initio parameters is the 
application of variational methods, such as the general- 
ized Gutzwillcr ansatz, which has recently been applied 
to the 3d-ferromagnet nickel^. 



VII. COMPARISON WITH OTHER AB-INITIO 
METHODS 

In this section we compare our ab-initio approach 
with existing and established first-principles methods for 
electronic-structure calculations . 

DFT-based approaches encounter the problem that the 
density functional of the exchange-correlation energy is 
not known, and hence additional assumptions and ansatz 
are necessary. Within LDA the assumption is made that 
the particle density to be determined is the same as that 
of an effective one-particle system, and the ansatz for the 
(local) density dependence of the exchange-correlation 
energy is chosen so that in the case of the homogeneous 
electron gas the most accurately known results are re- 
produced. So results of many-body theory for the ho- 
mogeneous electron gas enter, but many-body theory is 
not directly applied for the lattice system (an inhomoge- 
neous electron gas) . In LDA the (local) Hartree potential 
for the inhomogeneous electron system is explicitly taken 
into account when solving the Kohn-Sham equations, but 
the exchange (Fock) contribution is only approximately 
considered (in a form which is correct only for the ho- 
mogeneous electron system). Therefore, whereas self- 
interaction terms exactly cancel, if the Hartree and Fock 
contributions are treated on the same level, they usually 
do not cancel in LDA. Our approach, on the other hand, 
is free from the problem of self-interaction, as can im- 
mediately be seen from (|45f) . because the Hartree- and 
exchange terms are treated on the same level. Further- 
more, many-body theory is directly applied for the lat- 
tice electron system instead of only for the homogeneous 
electron gas. From the very beginning we never assume 
a constant electron density or a dependence on the local 
density. The effects which the "gradient correction" and 
nonlocal density schemes corrections to LDA aim at are 
automatically considered in our approach. 

The standard Hartree-Fock approximation in first 
quantization requires the solution of the equations 

(-|V + V(r)W) + £ / dV^^(r) 

jcr' 

^ f „ ,e 2 c9*(r'Wr') 
- E J rfV jr .Vj W = em(*) ■ 

For equal spin and i — j the second (Hartree) and the 
third (Fock or exchange) contributions exactly cancel 
each other so that no self-interaction occurs. The aim 
of the recent self-interaction corrected (exact exchange) 



schemes^ as LDA improvements is, therefore, to consider 
these exchange terms rigorously. But the effective ex- 
change (Fock) potential is a non-local potential, and the 
one-particle wave-function ifi to be determined enters 
not only at the position r but also at all other positions 
r'. For that reason HFA treatments in first quantization 
and exact exchange schemes are computationally more 
complicated than LDA methods, especially for a peri- 
odic solid. Within our approach HFA- calculations are 
as simple as a Hartree calculation, since only expecta- 
tion values have to be evaluated. But in our method 
self-interaction terms also exactly cancel. Therefore, our 
approach should be as good as the exact-exchange ap- 
proach but easier in practical applications. 

When using an RPA (screened HFA) as the many- 
body approximation, our approach should be similar to 
the GWA. But in GWA one usually uses as the Green 
function (G-line) the LDA (effective one-particle) result 
whereas we suggest using a selfconsistently calculated 
Green function (within the many-body approximation 
used). Also, when using the LDA Green function, only 
a part of the exchange is taken into account unlike the 
exact exchange of our approximation. 

All the recent attempts to combine ab-initio and many- 
body methods (on different levels concerning the many- 
body treatment), namely, the LDA+U, the LDA++, and 
the LDA+DMFT-schemes, have in common that they 
start from an LDA ab-initio calculation and obtain their 
one-particle bands and density of states from this LDA 
treatment, and then add an interaction (correlation) term 
with an on-site Hubbard-U term. Then, on one hand, 
some of the correlations and interactions are already in- 
cluded in the effective one-particle energies, namely those 
on the LDA (homogeneous electron gas) level. On the 
other hand, since other correlations are treated in many- 
body theory for the Hubbard U, there may be a double 
counting of interaction (correlation) contributions, which 
is hard to justify and the magnitude of which is difficult 
to estimate. Furthermore, in these theories the Hubbard 
U usually is an additional free parameter, so that these 
are not really "ab-initio" (first-principles) approaches. In 
the our approach we do not start from the LDA but from 
an effective band structure calculation on the Hartree 
level. Therefore no correlation terms are implicitly in- 
cluded within the one-particle band structures and the 
problem of double counting of interaction terms does not 
occur. Furthermore, all Coulomb matrix elements are 
calculated from "first principles" . 



VIII. SUMMARY AND OUTLOOK 

In this paper we have suggested a new approach for 
combining ab-initio and many-body methods for the cal- 
culation of the electronic properties of solids. The start- 
ing point is a traditional band-structure calculation for 
an effective (auxiliary) one-particle Hamiltonian, which 
can be the Hartree-Hamiltonian. This yields, in par- 
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ticular, the eigenfunctions in the form of Bloch func- 
tions. Keeping only a finite number of J band indices 
restricts and truncates the Hilbert space for further cal- 
culations. We use the Marzari-Vanderbilt algorithm to 
construct maximally localized Wannier functions (within 
the truncated one-particle Hilbert space). Then all the 
one-particle (tight-binding) and two-particle (Coulomb) 
matrix elements between these Wannier functions can be 
calculated. The strong localization guarantees that only 
on-site matrix elements and near-neighbor inter-site ma- 
trix elements have to be calculated. We are left with a 
many-body Hamiltonian in second quantization but with 
parameters determined from first principles for any given 
material. This should be solved by a suitable many-body 
approximation. The simplest approximation is the (un- 
screened) HFA, but improved methods (e.g., RPA or non- 
perturbative many-body approximations like DMFT) are 
also feasible. Our suggested approach is free from the 
problems of double counting of correlation effects and 
self-interaction and considers exchange contributions ex- 



actly. It does not rely on assumptions based on the homo- 
geneous electron gas or a dependence on the local electron 
density density. An inhomogeneous (lattice) electron sys- 
tem is considered right from the beginning. 

An application of this scheme to the 3d transition met- 
als Fe, Co, Ni and Cu using the unscreened HFA as the 
many-body method is presented in the following paper. 34 



Acknowledgments 

This work has been supported by a grant from the 
Deutsche Forschungsgemeinschaft No. Cz/31-12-1. It 
was also partially supported by the Department of En- 
ergy under contract W-7405-ENG-36. This research 
used resources of the National Energy Research Scien- 
tific Computing Center, which is supported by the Office 
of Science of the U.S. Department of Energy under Con- 
tract No. DE-AC03-76SF00098. 



1 P. Hohenberg, W. Kohn, Phys. Rev. 136, B 864 (1964) 

2 R.M. Dreizler, E.K.U. Gross, Density Functional Theory, 
Springer (Berlin, Heidelberg, New York 1990) 

3 H. Eschrig, The Fundamentals of Density Functional The- 
ory (Teubner Stuttgart, Leipzig 1996) 

4 W. Kohn, L.J. Sham, Phys. Rev. 140, A1133 (1965) 

5 J.H. Rose, JR. Smith, F. Guinea and J. Ferrante, Phys. 
Rev. B 29, 2963 (1984). 

6 Chapter 7 of Ref. □ 

7 O. Gunnarsson, R.O. Jones, Physica Script. 21, 394 (1980) 

8 M. Stadele, J.A. Majewski, P.Vogl, A. Gorling, Phys. Rev. 
Letters 79, 2089 (1997) 

9 L. Hedin, B.I. Lundquist, Solid State Physics 23, p.l (Eds.: 
F. Seitz, D. Turnbull, H. Ehrenreich, Academic Press 1969) 

10 For a recent review see: W. Aulbur, L. Jonsson, J.W. 
Wilkins in: Solid State Physics 54, p. 2 (Eds.: H. Ehren- 
reich, F. Spaepa, Academic Press 2000) 

11 V.I. Anisimov, J. Zaanen, O.K.Andersen, Phys. Rev. B 44, 
943 (1991) 

12 M.M. Steiner, R.C. Albers, D.J. Scalapino, L.J. Sham, 
Phys. Rev. B 43, 1637 (1991) 

13 M.M. Steiner, R.C. Albers, L.J. Sham, Phys. Rev. B 45, 
13272 (1992); Phys. Rev. Lett. 72, 2923 (1994) 

14 V.I. Anisimov, A.I. Poteryaev, M.A. Korotin, A.O. 
Anokhin, G. Kotliar, J. Phys. Cond. Matter 9, 7359 (1997) 

15 A.I. Lichtenstein, M.I. Katsnelson, Phys. Rev. B 57, 6884 
(1998) 

16 V. Drchal, V. Janis, J. Kudrnovsky, Phys. Rev. B 60, 
15664 (1999) 

17 M.I. Katsnelson, A.I. Lichtenstein, J. Phys. Cond. Matter 
11, 1037 (1999) 

18 A. Liebsch, A. Lichtenstein, Phys. Rev. Lett. 84, 1591 
(2000) 

19 T. Wegner, M. Potthoff, W. Nolting, Phys. Rev. B 61, 
1386 (2000) 

20 LA. Nekrasov, K. Held, N. Bliimer, A.I. Poteryaev, V.I. 
Anisimov, D. Vollhardt, Eur. Phys. J. B 18, 55 (2000) 



21 For a review see: A. George, G. Kotliar, W. Krauth, M.J. 
Rozenberg, Rev. Mod. Phys. 68, 13 (1996) 

22 K. Held, LA. Nekrasov, G. Keller, V. Eyert, N. Bliimer, 
A.K. McMahan, R.T. Scalettar, T. Pruschke, V.I. Anisi- 
mov, D. Vollhardt, in: Quantum Simulations of Complex 
Many-Body systems: From Theory to Algorithms, p. 175, 
(eds.: J. Grotendorst, D. Marx, A. Muramatsu, NIC Series 
Vol. 10, Forschungszentrum Jiilich 2002) 

23 N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 
(1997). 

24 O.K. Andersen, Phys. Rev. B 12, 3060 (1975) 

25 H.L. Skriver, The LMTO Method (Springer- Verlag, Heidel- 
berg 1984) 

26 D. Singh, Plane waves, pseudopotentials and the LAPW 
method (Kluwer Academic, Amsterdam 1994) 

27 P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, J. 
Luitz, WIEN2k, An Augmented Plane Wave + Local Or- 
bitals Program for Calculating Crystal Properties (Karl- 
heinz Schwarz, Techn. Universitat Wien, Austria), 2001. 
ISBN 3-9501031-1-2 

28 K. Koepernik and H. Eschrig, Phys. Rev. B 59, 1743 
(1999); I. Opahle, K. Koepernik, and H. Eschrig, Phys. Rev. 
B 60, 14035 (1999). 

29 I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 
035109 (2002) 

30 W. Ku, H. Rosner, W. E. Pickett, and R. T. Scalettar, 
Phys. Rev. Lett. 89, 167204 (2002) 

31 I. Schnell, G. Czycholl, R.C. Albers, Phys. Rev. B 65, 
075103 (2002) 

32 GD. Mahan, Many-Particle Physics (Plenum Press New 
York 1990) 

33 W. Weber, J. Biinemann, F. Gebhard, inBand- 
Ferromagnetsim (Lecture Notes in Physics, Vol. 580, p. 9 
(eds.: K.Baberschke, M. Donath, W. Nolting, Springer 
Berlin 2001); J. Biinemann, W. Weber, F. Gebhard, Phys. 
Rev. B 57, 6896 (1998) 

34 I. Schnell, G. Czycholl, R.C. Albers, Phys. Rev. B (2003), 



11 

following paper 



